#include "CutCallback.h"
#include <math.h>

int CutCallback::numberOfAddedViolatedInequalities = 0;

CutCallback::CutCallback( IloEnv& _env, string _cut_type, double _eps, Instance& _instance, IloBoolVarArray& _x, IloBoolVarArray& _z, IloBoolVar** _e, unsigned _k  ) :
	LazyConsI( _env ), UserCutI( _env ), env( _env ), cut_type( _cut_type ), eps( _eps ), instance( _instance ), x( _x ), z( _z ), e( _e ), k( _k )
{
	arc_weights.resize( 2 * instance.n_edges );
	if (cut_type=="dcc")
	{
		list<pair<u_int,u_int>> arcs;
		for (Instance::Edge edge : instance.edges)
		{
			u_int v1 = min(edge.v1,edge.v2);
			u_int v2 = max(edge.v1,edge.v2);
			
			arcs.push_back(pair<u_int,u_int>(v1,v2));
			if (v1 != 0)
			{
				arcs.push_back(pair<u_int,u_int>(v2,v1));
			}
		}
		u_int n=instance.n_nodes;
		u_int m=instance.n_edges;
		m = 2*m-(n-1);
		maxFlow = new Maxflow(instance.n_nodes,m,arcs);
	}
}

CutCallback::~CutCallback()
{
	delete maxFlow;
}

void CutCallback::separate()
{
	if( cut_type == "dcc" ) connectionCuts();
	else if( cut_type == "cec" ) cycleEliminationCuts();
}

/*
 * separation of directed connection cut inequalities
 */
void CutCallback::connectionCuts()
{
	if (!lazy)
		return;

	//cout << "----- CALLED CUT CALLBACK ------" << endl;

	u_int n = instance.n_nodes;
	u_int m = instance.n_edges;

	IloNumArray zval(env,n);
	LazyConsI::getValues(zval,z);
	double* zvalCopy = new double[n];
	int* cut = new int[n];

	u_int numFractionalNodes = 0;
	for (int i=0;i<n;i++)
	{
		zvalCopy[i] = zval[i];
		if ((zval[i] > eps && zval[i] < 1-eps) || zval[i] > 1+eps)
			numFractionalNodes++;
	}
	//cout << "Num fractional nodes: " << numFractionalNodes << endl;

	u_int numFractionalEdges = 0;
	double* capacity = new double[2*m-(n-1)];
	try
	{
		u_int capIndex = 0;
		for (Instance::Edge edge : instance.edges)
		{
			u_int v1 = min(edge.v1,edge.v2);
			u_int v2 = max(edge.v1,edge.v2);
			
			double temp = (double) LazyConsI::getValue(*e[v1*n+v2]);
			if ((temp > eps && temp < 1-eps) || temp > 1+eps)
				numFractionalEdges++;

			capacity[capIndex++] = temp;
			if (v1 != 0)
			{
				temp = (double) LazyConsI::getValue(*e[v2*n+v1]);
				if ((temp > eps && temp < 1-eps) || temp > 1+eps)
				numFractionalEdges++;

				capacity[capIndex++] = temp;
			}
		}
		//cout << "Num fractional edges: " << numFractionalEdges << endl;

		maxFlow->update(0,1,capacity);

		u_int numCutValueSmallerNodeValue = 0;
		for (u_int i=1;i<n;i++)
		{
			double nodeValue = zval[i];
			maxFlow->update(0,i);
			double cutValue = maxFlow->min_cut(3, cut);

			if (cutValue < nodeValue)
			{
				//cout << "Found cut: " << cutValue << " < " << nodeValue << endl;
				numCutValueSmallerNodeValue++;
				u_int cutNodeTypes[3] = {0,0,0};
				for (u_int j=0;j<n;j++)
				{
					cutNodeTypes[cut[j]]++;
				}
				//for (u_int j=0;j<3;j++)
				//{
				//	cout << "  type " << j << ": " << cutNodeTypes[j] << endl;
				//}
				//if (cut[0] == 0)
				//	cout << "  Node 0 is type 0!" << endl;
				//if (cut[i] != 2)
				//	cout << "  Node i=" << i << " is type " << cut[i] << endl;

				IloExpr minCutExpr(env);

				for (Instance::Edge edge : instance.edges)
				{
					u_int v1 = min(edge.v1,edge.v2);
					u_int v2 = max(edge.v1,edge.v2);
			
					if (cut[v1] != 2 && cut[v2] == 2)
						minCutExpr += *e[v1*n+v2];
					if (v1 != 0)
					{
						if (cut[v2] != 2 && cut[v1] == 2)
							minCutExpr += *e[v2*n+v1];
					}
				}

				LazyConsI::add(minCutExpr >= z[i]);
				numberOfAddedViolatedInequalities++;
			}
		}
		//cout << "Num cut value smaller node value: " << numCutValueSmallerNodeValue << endl;
	}
	catch( IloException& e ) {
		cerr << "CutCallback: exception " << e.getMessage();
		exit( -1 );
	}
	catch( ... ) {
		cerr << "CutCallback: unknown exception.\n";
		exit( -1 );
	}

	//cout << "------- END CUT CALLBACK ------------" << endl << endl;
}

/*
 * separation of cycle elimination cut inequalities
 */
void CutCallback::cycleEliminationCuts()
{
	if (!lazy)
		return;

	try {

		u_int n = instance.n_nodes;
		u_int m = instance.n_edges;

		//IloNumArray xval( env, 2 * m );
		//IloNumArray zval( env, n );

		vector<int> sources, targets, arcIndices;

		for (unsigned ei=0; ei<m; ei++)
		{
			Instance::Edge& edge = instance.edges[ei];

			unsigned vi[2] = {edge.v1, edge.v2};
			for (unsigned i=0;i<2;i++)
			{
				unsigned v1 = vi[0^i];
				unsigned v2 = vi[1^i];

				if (v2 != 0 && v1 != 0)
				{
					IloBoolVar* edgeVar = e[v1*n+v2];
					IloNum ret;
					if( lazy )
					{
						ret = LazyConsI::getValue(*edgeVar);
					}
					else
					{
						ret = UserCutI::getValue(*edgeVar);
					}
					//cout << ret << endl;
					double newWeight = 1.0-ret;
					arc_weights[ei+i*m] = newWeight;

					if (newWeight <= 1-eps)
					{
						sources.push_back(v1);
						targets.push_back(v2);
						arcIndices.push_back(ei+i*m);
					}
				}
				else
				{
					arc_weights[ei+i*m] = 2*k;
				}
			}
		}

		//cout << "Adding cycle: ";
		for (unsigned i=0;i<sources.size();i++)
		{
			SPResultT t = shortestPath(targets[i], sources[i]);
			u_int numEdgesInCycle = t.path.size();
			double myWeight = arc_weights[arcIndices[i]];
			if (t.weight + myWeight < 0.8 - numEdgesInCycle*eps)
			{
				IloExpr sumCycle(env);
				IloBoolVar* firstArcVar = e[sources[i]*n+targets[i]];
				//cout << firstArcVar->getName() << ", ";
				sumCycle += *firstArcVar;
				for (unsigned arcIndex : t.path)
				{
					int swap = (arcIndex >= m);
					int edgeIndex = arcIndex - m*swap;
					Instance::Edge& edge = instance.edges[edgeIndex];
					unsigned vi[2] = {edge.v1, edge.v2};
					IloBoolVar* arcVar = e[vi[0^swap]*n+vi[1^swap]];
					//cout << arcVar->getName() << ", ";
					sumCycle += *arcVar;
				}
				//cout << endl;

				IloNum pathSize = (IloNum)t.path.size();
				if( lazy )
					LazyConsI::add( sumCycle <= pathSize + pathSize*eps);
				else
					UserCutI::add( sumCycle <= pathSize + pathSize*eps);

				numberOfAddedViolatedInequalities++;
				sumCycle.end();

			break;
			}
		}

		/*if( lazy ) {
			LazyConsI::getValues( xval, x );
			LazyConsI::getValues( zval, z );
		}
		else {
			UserCutI::getValues( xval, x );
			UserCutI::getValues( zval, z );
		}*/



		// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
		// TODO find violated directed connection cut inequalities
		// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++

		// add found violated cut to model
		//if( lazy ) LazyConsI::add( ... );
		//else UserCutI::add( ... );

		//xval.end();
		//zval.end();

	}
	catch( IloException& e ) {
		cerr << "CutCallback: exception " << e.getMessage();
		exit( -1 );
	}
	catch( ... ) {
		cerr << "CutCallback: unknown exception.\n";
		exit( -1 );
	}
}

/*
 * Dijkstra's algorithm to find a shortest path
 * (slow implementation in time O(n^2))
 */
CutCallback::SPResultT CutCallback::shortestPath( u_int source, u_int target )
{
	u_int n = instance.n_nodes;
	u_int m = instance.n_edges;
	vector<SPNodeT> nodes( n );
	vector<bool> finished( n, false ); // indicates finished nodes

	// initialization
	for( u_int v = 0; v < n; v++ ) {
		nodes[v].pred = -1;
		nodes[v].pred_arc_id = -1;
		if( v == source ) nodes[v].weight = 0;
		else nodes[v].weight = numeric_limits<double>::max();
	}

	while( true ) {

		// find unfinished node with minimum weight to examine next
		// (should usually be done with heap or similar data structures)
		double wmin = numeric_limits<double>::max();
		u_int v;
		for( u_int u = 0; u < n; u++ ) {
			if( !finished[u] && nodes[u].weight < wmin ) {
				wmin = nodes[u].weight;
				v = u;
			}
		}

		// if all reachable nodes are finished
		// or target node is reached -> stop
		if( wmin == numeric_limits<double>::max() || v == target ) break;

		// this node is finished now
		finished[v] = true;

		// update all adjacent nodes on outgoing arcs
		list<u_int>::iterator it;
		for( it = instance.incidentEdges[v].begin(); it != instance.incidentEdges[v].end(); it++ ) {
			u_int e = *it; // edge id
			u_int a; // according arc id
			u_int u; // adjacent node
			if( instance.edges[e].v1 == v ) {
				a = e;
				u = instance.edges[e].v2;
			}
			else {
				a = e + m;
				u = instance.edges[e].v1;
			}
			// only examine adjacent node if unfinished
			if( !finished[u] ) {
				// check if weight at node u can be decreased
				if( nodes[u].weight > nodes[v].weight + arc_weights[a] ) {
					nodes[u].weight = nodes[v].weight + arc_weights[a];
					nodes[u].pred = v;
					nodes[u].pred_arc_id = a;
				}
			}
		}
	}

	SPResultT sp;
	sp.weight = 0;
	int v = target;
	while( v != (int) source && v != -1 ) {
		int a = nodes[v].pred_arc_id;
		if( a < 0 ) break;
		sp.weight += arc_weights[a];
		sp.path.push_back( a );
		v = nodes[v].pred;
	}
	return sp;
}

